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Abstract: A program to simulate the production of heavy quarks through the 
boson-gluon fusion process in e^p collisions is presented. The full electroweak 
structure of the electron-gluon interaction is taken into account as well as the 
masses of the produced heavy quarks. Higher order QCD radiation is treated 
using initial and final state parton showers, and hadronization is performed 
using the Lund string model. Physics and programming aspects are described 
in this manual. 



1 Physics of included processes 

The production of heavy quarks in ep collisions is dominated by boson-gluon fusion (BGF) 
into a heavy quark-antiquark pair 

V(q) + g(p)^Qf(pf) + Qr(p f >) , (l) 

occuring as 0(a s a 2 ) parton level subprocess in the ep scattering process 

e ± (l e ) + p(P) - l\V) + Qfipf) + Q f (p f ) + X . (2) 

Here the symbols in brackets denote the corresponding four-momenta. In the charged 
current (CC) case V is the IV ± -boson whereas in the neutral current (NC) case it corre- 
sponds to 7/Z exchange and the produced quark and antiquark have the same flavour 
/. Beside the normal DIS variables defined by 

Q 2 = -q 2 = ~{h - 1') 2 , W^(P + qf , * = , V = J^> ( 3 ) 

out of which only two are independent, three additional independent variables are needed 
to completely specify the process in eq. @. These are here taken as the momentum 
fraction x g of the gluon relative to the proton, i.e. p = x g P, the variable 



P-q 



(4) 



1 Heisenberg fellow 
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and the azimuthal (around the boson axis) angle $ between the lepton and hadron planes 

(p x C) • (p x p}) 

008$ = =r — — (5) 

| p x Z e 1 1 p x p") i 

measured in the boson-gluon CM frame, i.e. p/ +p/' = 0. The gluon momentum fraction 
is related to the usual Bjorken-x variable by 

s 

x„ = x H > x . (6) 

ys 

The variable 2 is related to the angle of the QQ'-a.xis with respect to the boson-gluon 
axis in this subsystem cms and gives the heavy quark transverse momentum through the 
relation p\ = sz(l — z) + z(mj, — m 2 ) — rrxh. As usual, s is taken to denote the invariant 
mass square of the QQ'-sribsystem, i. e . s = (p/ +p/') 2 - 

The cross section for process (0) is then given by 

a{e ± p -> QQ'X) = J dy J dQ 2 J dx g J dz f d<$> g(x g , M 2 g ) h(y, Q 2 , x g , z, $) (7) 

which is a convolution of the gluon density g(x g , M 2 ) and a QCD part h for the subprocess. 
The latter has been calculated in ref. taking proper account of the heavy quark masses 
and the complete electroweak structure for both charged and neutral current processes 
including the 7 — Z° interference and allowing for arbitrary longitudinally polarization of 
the e beam. Both the deep inelastic and photoproduction (Q 2 — > 0) region is covered. 

For the numerical evaluation [| of this cross section the electroweak parameters are 
taken from Lepto 6.5 The Kobayashi-Maskawa matrix elements have the default val- 
ues of Lepto 6.5 and the quark masses the default values of Jetset 7.4 0], in particular 
m c = 1.35 GeV and nib = 5 GeV. All these default values can be changed by the user. The 
mass scales for the gluon density and for a S) M g and M s , are taken to be M g = M s = 
the number of flavours, Nf, and Aqqd used in the first order expression for a s are taken 
from Jetset which means that Aqcd is taken from the parton densities. A selection 
of parametrizations for the gluon density in the proton are provided as alternatives (cf. 
LST(15) 0). The numerical uncertainties in the cross sections are discussed in [Q. 

The next-to-leading order (NLO) QCD corrections to the inclusive heavy quark dis- 
tributions (rapidity and p±) have been calculated for photoproduction (Q 2 <C A 2 ) || and 
recently also for deep- inelastic scattering (DIS), i.e. Q 2 3> A 2 M. In the case of pho- 
toproduction, the 0(a 2 ) corrected distributions are found to be similar in shape to the 
leading-order ones, i.e. they can be described by an approximately constant i^-factor. It 
is therefore expected that the distributions from Aroma describe the shapes of inclusive 
heavy quark distributions reasonably well, whereas the absolute normalization is under- 
estimated [Q. In contrast, the ii'-factor of DIS is strongly kinematics dependent and 
therefore not easily taken into account. In addition, the available NLO result only applies 
for inclusive quantities and not for the exclusive final state needed for MC simulation. 

For a proper description of the complete event topology higher order corrections are 
important. Gluon bremsstrahlung off the heavy quark and antiquark can be approxi- 
mately taken into account by the use of parton shower (PS) simulation algorithms. From 
the experience of jets in e + e~ annihilation one expects this approach, which includes the 
leading logarithm contribution of higher corrections, to provide a better description of 
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detailed jet properties, such as hardness and width, as compared to exact next-to-leading 
order matrix elements. Such a shower algorithm [|J is therefore applied, but one should 
note that the rate of clearly separable additional jets need not be quite correct. 

Production of heavy quarks in lowest order, eq. (1), is characterized by a pair of heavy 
quarks that are essentially back-to-back in the boson-proton CMS and with small overall 
transverse momentum p±(QQ). Small deviations from this can arise from a non- vanishing 
primordial k± of the incoming gluon, which is here generated according to a Gaussian 
distribution. Larger deviations from this topology may arise through parton radiation 
from the gluon entering the fusion process. Such initial state radiation also affects the 
inclusive heavy quark spectra, which in NLO accuracy follow rather closely those at the 
Born level (at least at low Q 2 ). This implies constraints on the initial radiation and 
therefore one does not expect large effects. Initial state parton showers are implemented 
by suitably modifying the algorithm of ref . |7| , in particular allowing for a smooth Q 2 — > 
limit. 

In ep scattering heavy quarks can also be produced through (i) mixing in the leading 
order charged current process e + q — > v + Q [EJ, (ii) the scattering on intrinsic charm 
quarks in the proton ||, (Hi) resolved photon contributions leading to hadronic production 
processes of heavy flavours (e.g. gg — > QQ) and (iv) diffractive heavy flavour excitation 
0. These processes, which are not included in Aroma, are expected to have much smaller 
cross sections than the BGF process in the HERA energy range. 

The Aroma program is also applicable to the case of light flavours. However, without 
the heavy quark mass as a cutoff against the divergences in the matrix element it is then 
necessary to apply a lower cut on the p± of the produced quarks in the boson-gluon fusion 
cms. The program thus allows to set a pj_-cutoff, which may also be applied in heavy 
quark generation. 



2 The Monte Carlo implementation 

The Monte Carlo simulation model is based on the following main ingredients: (i) The 
complete matrix elements to order a 2 a s for the boson-gluon fusion process, eq. (0), 
(including the masses of the produced heavy quarks and the full electroweak structure of 
the interaction), (ii) gluon emission from the QQ' system in a parton shower approach 
HI, (Hi) initial state parton shower, (iv) possible soft colour interactions as a mechanism 
for rapidity gaps [IIJ, (v) hadronization with the Lund string model [II], |J and heavy 
flavour decay. 



The importance sampling method in Divonne |12| is used to generate phase space 
points according to the differential cross-section formula in eq. (|7|). From the phase space 
point (y, Q 2 , x g , z, $), the four-momenta of all partons (particles) are calculated. 

Similarly to the qq pair produced in e + e~ annihilation the heavy QQ' pair can emit 
bremsstrahlung gluons (some of which may split perturbatively into gluon or qq pairs) 
thereby creating a cascade, or shower, at the parton level. To simulate this the model 
for e + e~ [|J is applied. In doing so an uncertainty arises as to what scale to use for 
the maximum off-shellness of the heavy quark and antiquark that initiate the shower. 
Different options are provided, see LSTHFL(7), with the default chosen as {in±Q + m ± Q/) 2 
since this measures the momentum transfer in similarity with the case in the gg — ► 
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QQ process in hadron collisions. To account for parton radiation from the incoming 
gluon, the initial state parton shower in Lepto has been incorporated using the scale 
min(/z 2 , W 2 - s, s + Q 2 /2) where u 2 is set by LSTHFL(9) (default is (m± Q + m L Q,) 2 ). 

The multiparton final state is then supplemented by the proton remnant. In case of an 
incoming gluon, the remnant contains the three valence quarks. The energy-momentum 
of the remnant is divided and a small relative transverse momentum introduced when 
this system is split |T3| , |3| into a quark and a diquark. Together with the produced heavy 
antiquark and quark, respectively, these spectator partons form two separate colour triplet 
strings. Before hadronising these strings with Jetset 7.4 |4| extra soft colour interactions 
may be taken into account^ which is regulated by LST(34). This results in a Monte Carlo 
model for heavy flavour production in ep collisions that generates complete events with 
the full information on both the parton and the hadron level. 



3 Description of program components 

The program is written in FORTRAN 77 and consists of a set of subroutines that form 
an add-on package to Lepto 6.5 B. In addition, Jetset 7.4 and Pythia 5.7 M and 



the Divonne program [0, which are all available via Cernlib, are used. The user has 
to ensure that the Aroma routines are loaded such that they replace some routines in 
Lepto with their modified versions, i.e. normally Aroma should be loaded first. All 
subroutine and common-block names start with HF to avoid name clashes. Exceptions 
are: AROMA and AROINI which the user calls; DFUN and DVNOPT which replaces 
subroutines in Lepto 6.5 and the block data ARODAT. 

The heavy flavour generation is made by calling AROMA. The parameters, cuts and 
switches for heavy flavour generation are in arrays PARHFL, CUTHFL and LSTHFL of 
common HFDATA described below. 



3.1 Subroutines and functions 



The following routines are called by the user: 
SUBROUTINE AROINI(LFILE,LEPIN,PLZ,PPZ,INTER) 

to initialize the event generation procedure, see also 



Purpose: 
Arguments: 
LFILE : 
LEPIN : 
PLZ, PPZ : 



INTER : 



= normally when using AROMA. 

type of lepton in Jetset 7.4 code, i.e. 11 = e~, —11 = e + . 
momentum (GeV/c) for incoming lepton and nucleon, respectively, along 
the z-axis (if both non-zero, i.e. colliding beams, they must have opposite 
signs) . 

type of interaction to be simulated, 
electromagnetic (EM), i.e. 7 exchange, 
weak charged current (CC), i.e. W ± exchange, 
weak neutral current, i.e. Z° exchange. 



2 This will give some production of qq resonances for sufficiently small masses of the qq system. How- 
ever, the treatment is rather crude and only J/t/j and T are produced for cc and bb respectively. 
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=4: neutral current (NC), i.e. j/Z exchange. 
SUBROUTINE AROMA 

Purpose: to administer the generation of one event of the kind specified by the last 
AROINI call. 

Remark: If the error flag LST(21) is non-zero then the event should be rejected 
(should only occur rarely). 

SUBROUTINE LFRAME(IFRAME,IPHI) 

Purpose: to transform the event between different frames, see ||. 

In the following list all subroutines (S) and functions (F) are briefly described. The 
order is as they appear in the code and reflects the flow in the program. Further details 
about their purpose and procedures used are given by comments in the code. 

Routine Purpose 

ARODAT Block data to give default values to all switches and parameters. 
AROINI (S) To initialize program, see above. 

HFINIT (S) Initiate constants, calculate integrated cross section, set-up grid for event 

generation. Called by AROINI. 
AROMA (S) To generate an event, see above. 

HFLGEN (S) Administer heavy flavour generation based on exact matrix elements, called 
by AROMA. 

DVNOPT (S) Substitutes block data for Divonne. 
DFUN (F) Integrand for Divonne, calls HFFUNC. 
HFFUNC (F) Differential cross section summed over flavours. 



z. 



HFVECT (S) Transforms XX(1), . . ., XX(4) to y, Q 2 , x 
HFCHCK (S) Check if y, Q 2 , x g and z are within bounds for given flavours. 
HFTOTD (F) Differential cross section for given flavours. 
HFVCT1 (S) Variable mapping for 7-exchange, called by HFVECT. 
HFVCT2 (S) Variable mapping for W- and Z-exchange, called by HFVECT. 
HFMEPS (S) Modified version of LMEPS in Lepto. 
HFSCAL (S) Modified version of LSCALE in Lepto. 
HFSSPA (S) Modified version of LYSSPA in Lepto. 
HFANGL (F) Modified version of ULANGL in Jetset. 
HFDBRB (S) Modified version of LUROBO in Jetset. 

The following subroutines in Divonne 4 are used: PARTN, INTGRL, RANGEN. 



3.2 Common blocks 

The common blocks intended for communication with the program are HFDATA and 
LEPTOU. Common LUJETS in Jetset 7.4 j| is used to store the event record and 
is therefore essential. All variables are given sensible default values in the block data 
ARODAT, as shown by (D=. . . ) below. Overwriting the default values should be made 
before calling AROINI. 
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COMMON /HFDATA/ PARHFL(10),CUTHFL(8),LSTHFL(10) 



Purpose: 

PARHFL(l) : 
PARHFL(2) : 
PARHFL(3) : 

PARHFL(4) : 
PARHFL(5) : 

PARHFL(6) : 



PARHFL(7-10) 



CUTHFL(l) 
CUTHFL(2) 
CUTHFL(3) 



CUTHFL(4) 
CUTHFL(5) 
CUTHFL(6) 
CUTHFL(7) 
CUTHFL(8) 

Remark: 



contains input switches and input parameters to specify physics, kine- 
matic cuts and numerical procedures, as well as output flags. Note that 
PARHFL and CUTHFL are in double precision. 

total cross section in pb for heavy flavour production. Calculated in the 
initialization taking the hard cuts in CUTHFL into account. 
(D=0.2) precision parameter SPRDMX for call to PARTN in the Dl- 
VONNE program; smaller values give finer partitioning. 
(D=1000.) parameter NPT for call to INTGRL in Divonne. Specifies 
maximum function evaluations in each subregion for integrating total 
heavy flavour cross section. 

(D=0. GeV 2 ) defines the minimum p\ of the quarks in the boson-gluon 
cms, must be used when producing light quarks. 

fraction of phase space points generated within the 'hard' cuts given in 
CUTHFL (see below) that are accepted by the 'soft' cuts in CUT (see 
common LEPTOU in ||) when generating the complete events, 
cross-section in pb for the generated event sample. Given by the product 
of PARHFL (1) and PARHFL (5), thus taking both hard and soft cuts 
into account. To be used at the end of event generation (updated for 
each event) to get absolute cross section normalisation of generated 
event sample, 
not used 

y mhi (D=0.) 

Umax (D = l.) 

Qmm (D=0.) (GeV 2 ), can be set to zero when only heavy flavour pro- 
duction is switched on (cf. LSTHFL(l)). 
Qmax (D=1.D9) (GeV 2 ) 



X 



g,min 



(D=0.) 



x 9 , max (D=0.9999999) 

Zmin (D=0.) 
^max (D 1.) 

These 'hard' cuts are on the basic variables in the heavy flavour matrix 
element and can therefore be efficiently applied to limit the region of 
phase space. (Kinematic phase space limits are explicitly calculated 
and applied if stronger than these cuts.) In contrast, the 'soft' cuts 
in CUT (common LEPTOU in ||) are used to reject chosen points or 
partly generated events, resulting in less efficiently applied cuts. See 
PARHFL (5), PARHFL(6) above. 



LSTHFL(l) : (D=2) regulates heavy flavour generation. 

=0: DIS events, no heavy flavour production. 
= 1: DIS events and heavy flavour production. 
=2: only heavy flavour production. 
Remark: When 'normal' LEPTO events and heavy flavour events are mixed in 
the same run, the total cross section is then the sum of the normal one 
(PARL(23) or PARL(24)) and the heavy flavour one (PARHFL(6)). 
To avoid double counting precautions have been taken in AROINI by 
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LSTHFL(2) 
LSTHFL(4) 



LSTHFL(5) 
LSTHFL(6) 



LSTHFL(7) 



LSTHFL(8) 



LSTHFL(9) 



setting switches such that the same flavours should not be produced 
in other processes; i.e. boson-gluon fusion process with massless matrix 
elements in Lepto (LST(13)), through parton showers (MSTJ(45), 
IPY(8)) and via sea quark density parametrizations (LST(12)). 

LSTHFL(3) : (D=4,5) lightest and heaviest flavour to be produced in NC 
heavy flavour production, can be set between 1 and 6 for d, u, s, c, b, t. 
: (D=4) specifies process simulated 

=1: pure 7 contribution. 

=2: W exchange. 

=3: pure Z contribution. 

=4: full NC contribution (7 + 7/Z + Z). 
: internal, > for e + , < for e~ . 

: (D=ll) regulates which flavour combinations to be included in CC 
interactions. Should be set to I cs + I01 c b + 100J is + 1000/^, where I cs 
etc. is to exclude (1 to include) cs (or cs) etc. 

: (D=4) regulates the scale in the final state parton showers in QQ sys- 
tem. Value is multiplied by PYPAR(25). 

Q\ 
W 2 , 

s, i.e. same treatment as for gg-event in Lepto 
(m ±Q + m ± g) 2 , 



171 ±q + 171 ±0 ' wnere m ±Q is the transverse mass of the heavy quark. 
(D=2) simulation of QCD effects in hadronic final state, see [[|. 
no parton showers 

both initial and final state parton showers 
initial state parton shower 
final state parton shower 

(D=4) choice of scale, i.e. maximum parton offshellness, in QCD initial 
state cascade. Value is multiplied by PYPAR(26). 

Q 2 , 



m 2 ! r> + m 2 . 



MQ ' ••"l_Qi 

Same treatment as for gg-event in Lepto. 
Remark: For the initial state shower the scale used is min(/z 2 , W 2 — s, s + Q 2 /2) 
where a 2 is the scale chosen above (Q 2 should not be used for photo- 
production). 

COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U 

Purpose: contains input switches and input parameters to specify physics, soft kine- 
matic cuts and numerical procedures, as well as output flags and output 
variables. 

Remark: only changes for Aroma are listed here, otherwise see 0. 
LST(8) : set by AROINI corresponding to LSTHFL(8). 
LST(17)=1: not usable with Aroma. 

LST(24)=5: current event is heavy quark event from AROMA. 

LST(35) : there is no special sea quark treatment in the initial cascade in Aroma. 



s 
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COMMON /HFLVAR/ SHAT,ETA,Z,PHIQQ,PT2Q,TM2Q1,TM2Q2,ZJ,IFLQ,IFLQB, 
& PHIFL(0:9,0:2),DBETAZ 

Purpose: Variables at the parton level, values obtained from simulation of cross 

section. Note that DBETAZ is of type double precision. 
SHAT,ETA,Z,PHIQQ : variables s, x g , z and $, see section 1. 
PT2Q : p\ for produced Q (and Q) in the boson-gluon cms, see section 1. 
TM2Q1,TM2Q2 : m\ for produced Q and Q in the boson-gluon cms. 

COMMON /LFLMIX/ CABIBO(4,4) 

Purpose: Contains the Cabbibo-Kobayashi-Maskawa matrix elements squared for 
flavour mixing, see [[|. 

COMMON /LUDAT2/ ...PMAS(500,4)... 
Purpose: Various parameters in Jetset 

PMAS(IFL,1) : mass for quark IFL=l-6 for d, u, s, c, b, t used in the BGF matrix element 

COMMON /LUJETS/ N,K(4000,5),P(4000,5),V(4000,5) 
Purpose: Contains the record of the generated event, see (4]]. 

Remark: The event record contains the full development of the event including a 
few extra, but inactive lines used for internal testing purposes. 

The following common-blocks are mainly for internal use: HFBL1, HFBL2, HFJETS, 
HFJRDM, HFMACH, HFSNCO. 

Divonne common-blocks used: BNDLMT, SAMPLE, PRINT, QUADRE. 

4 Usage and Availability 

Aroma 2.2 should be loaded together with Divonne 0, Lepto 6.5 [[§, Jetset 7.4 
and Pythia 5.7 The program is a slave system, which the user must call from his 
own steering program. Information about the program, the source code and a demonstra- 
tion job can be obtained from the authors or via the WWW on the Aroma homepage, 
http : //www3 . tsl .uu. se/thep/aroma/. 
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